import pandas as pd
import numpy as np
import seaborn as sns
from sklearn.cluster import KMeans,AgglomerativeClustering,OPTICS
from sklearn.metrics import silhouette_score,davies_bouldin_score
from matplotlib import pyplot as plt
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
import matplotlib.cm as cm
data = pd.read_csv("data_parsed.csv")
data = data[data.columns[1:]]
tsnedata = TSNE().fit_transform(data)
_scores = [kk.score(data) * -1 for kk in
[k.fit(data) for k in
[KMeans(n_clusters=n,random_state=0) for n in range(8,100,10)]
]
]
x_ticks = list(range(8,100,10))
plt.plot(x_ticks, _scores, 'x-')
plt.xlabel('k')
plt.ylabel('Within-cluster sum of squares')
plt.title('The wcss score for each k')
plt.show()
Przybliżamy do dokładniejszego zbadania region 25-50, ponieważ jest w nim największe załamanie.
kmeansmap = {n: KMeans(n_clusters=n,random_state=0) for n in range(25,50,1)}
for val in kmeansmap.values():
val.fit(data)
_scores = [-val.score(data) for val in kmeansmap.values()]
x_ticks = list(range(25,50,1))
plt.plot(x_ticks, _scores, 'x-')
plt.xlabel('k')
plt.ylabel('Within-cluster sum of squares')
plt.title('The wcss score for each k')
plt.show()
x = kmeansmap.keys()
y = [silhouette_score(data,k.predict(data)) for (n,k) in kmeansmap.items()]
plt.figure(figsize=(16,8))
plt.xlabel('k')
plt.ylabel('score')
plt.title('The silhouette score for each k')
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f6514a365c0>]
y = [-davies_bouldin_score(data,k.predict(data)) for (n,k) in kmeansmap.items()]
plt.figure(figsize=(16,8))
plt.xlabel('k')
plt.ylabel('score')
plt.title('The Davies_Bouldin score for each k')
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f650c847978>]
Na podstawie tych danych, stwierdziliśmy, że k=45 jest najlepszą liczbą zgrupowań dla tych danych.
agglomemap = {n: AgglomerativeClustering(n_clusters=n) for n in range(25,50,1)}
for val in agglomemap.values():
val.fit(data)
x = kmeansmap.keys()
y = [silhouette_score(data,k.labels_) for (n,k) in agglomemap.items()]
plt.figure(figsize=(16,8))
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f650c8d2c18>]
y = [-davies_bouldin_score(data,k.labels_) for (n,k) in agglomemap.items()]
plt.figure(figsize=(16,8))
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f650c8c09b0>]
agglomemap2 = {n: AgglomerativeClustering(n_clusters=n) for n in range(46,56,1)}
for val in agglomemap2.values():
val.fit(data)
Sprawdzamy, czy n=47 nie jest tylko lokalnym minimum, poprzedzającym znacząco lepsze minimum, obliczając score dla n w zakresie 46-56
x = agglomemap2.keys()
y = [silhouette_score(data,k.labels_) for (n,k) in agglomemap2.items()]
plt.figure(figsize=(16,8))
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f6514a51358>]
y = [-davies_bouldin_score(data,k.labels_) for (n,k) in agglomemap2.items()]
plt.figure(figsize=(16,8))
plt.plot(x,y,'x-')
[<matplotlib.lines.Line2D at 0x7f650c857cf8>]
W obu przypadkach 47 to n dla którego zostaje osiągnięty najlepszy wynik
o = OPTICS().fit(data)
print(silhouette_score(data,o.labels_))
print(-davies_bouldin_score(data,o.labels_))
/home/mtizim/.local/lib/python3.6/site-packages/sklearn/cluster/_optics.py:803: RuntimeWarning: divide by zero encountered in true_divide ratio = reachability_plot[:-1] / reachability_plot[1:] -0.48469136652664857 -1.1860217484077489
OPTICS osiąga słabe wyniki na poprzednich testach, ale nawet w przypadku osiągnięcia lepszych wyników należałoby go odrzucić, co widać na poniższych wykresach porównujących otrzymane zgrupowania:
plt.scatter(tsnedata[:,0], tsnedata[:,1], c=kmeansmap[45].predict(data))
<matplotlib.collections.PathCollection at 0x7f650c357f60>
plt.scatter(tsnedata[:,0], tsnedata[:,1], c=agglomemap[47].labels_)
<matplotlib.collections.PathCollection at 0x7f650c037f98>
plt.scatter(tsnedata[:,0], tsnedata[:,1], c=o.labels_)
<matplotlib.collections.PathCollection at 0x7f650c2b62b0>
Bardzo dużo punktów nie zostało zgrupowanych przez OPTICS, przez co nie jest on dobrym algorytmem w tym przypadku
print(f"Niezgrupowane rekordy: {sum( 1 for i in o.labels_ if i == -1)}/{len(o.labels_)} ({sum( 100 for i in o.labels_ if i == -1)/len(o.labels_):.4f}%)")
Niezgrupowane rekordy: 10189/12330 (82.6358%)
Poprzednie wizualizacje pokazywały interesujące rezultaty na wykresie PCA, więc postanowiliśmy to porównać z naszymi modelami.
pca = PCA(n_components=2).fit_transform(data)
plt.scatter(pca[:,0], pca[:,1],c=agglomemap[47].labels_)
plt.show()
pca = PCA(n_components=2).fit_transform(data)
plt.scatter(pca[:,0], pca[:,1],c=kmeansmap[45].predict(data))
plt.show()
Widać, że widoczne 10 'klas' otrzymanych z PCA nie oddaje pełnego obrazu sytuacji.
pca = PCA().fit(data)
plt.figure(figsize=(9,6))
plt.plot(range(1, len(pca.explained_variance_ratio_)+1), np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
Żeby PCA sensownie obrazowało dane, musiałoby mieć przynajmniej 40 wymiarów, co nie jest wizualizowalne.